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Despite their topological complexity almost all functional properties of metabolic networks can 
be derived from steady-state dynamics. Indeed, many theoretical investigations (like flux-balance 
analysis) rely on extracting function from steady states. This leads to the interesting question, how 
metabolic networks avoid complex dynamics and maintain a steady-state behavior. Here, we expose 
metabolic network topologies to binary dynamics generated by simple local rules. We find that the 
networks' response is highly specific: Complex dynamics are systematically reduced on metabolic 
networks compared to randomized networks with identical degree sequences. Already small topolog- 
ical modifications substantially enhance the capacity of a network to host complex dynamic behavior 
and thus reduce its regularizing potential. This exceptionally pronounced regularization of dynamics 
encoded in the topology may explain, why steady-state behavior is ubiquitous in metabolism. 
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I. INTRODUCTION 

The general notion of network biology [l| proposes an 
abstract view on biological systems: The components' 
complex interaction pattern is represented by a mathe- 
matical graph consisting of nodes and links. The aim is 
to understand universal features and design principles 
of complex biological networks at a system-wide level 
0, [H, [|| . Recent findings include the ubiquity of heavy- 
tail degree distributions a 'bow-tie' structure of cer- 
tain network types 0, [f| , the presence of modules [1, 0, [I| 
and a similarity in motif content of functionally related 
networks @ . An important focus of research is to develop 
models of graph construction, whichyield similar statis- 
tical properties as the real graphs 0, El, E3| ■ This 
modeling task is complicated by the fact that dynamic 
performance is a criterion in the evolutionary shaping of 
some types of biological networks [TU, [l4| • Consequently, 
a current challenge is to incorporate dynamics into this 
general framework, i.e. to link topology and dynamic 
function. For gene regulatory networks huge progress 
has been made in the last few years in that regard: The 
motif content of eukaryotic genetic networks [9( has been 
shown [l3l | to correlate with the dynamic robustness pro- 
file of these motifs; the dominant branch of the largest 
attractor of a reduced yeast cell cycle network under the 
framework of Boolean dynamics coincides with the ex- 
perimentally observed sequence of cell states PH . 
In the case of metabolic networks the situation is far more 
involved. Since the dynamics of the metabolic concen- 
trations lack an approximately binary behavior, refined 
kinetic models have been developed [see, e.g., [HI, Il7j |. 
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Due to the limited knowledge of kinetic parameters, such 
models are often restricted to sub-networks, lacking the 
large-scale perspective, unless strong simplifications and 
abstractions are introduced in the dynamics. Many ap- 
proaches exploit an intriguing feature of metabolic net- 
work dynamics: the convergence to a steady state. Flux 
balance analysis (FBA) [lg, EH, for instance, retains 
the stoichiometric matrix (and, therefore, the interac- 
tion pattern of the network) to formulate hypotheses 
on the overall performance of the system, together with 
an objective function incorporating constraints on the 
metabolic dynamics. This method and related steady- 
state approaches have been successfully applied to the 
phenotypic prediction of various wild-type species and 
mutants under different environmental conditions [see, 

e.g., m m m. 

Here, we address a fundamental topological question: 
How well is a particular network designed to suppress 
complex dynamics? We select rules, which lead to 
transient (i.e., non steady-state) binary dynamics on 
metabolic networks and compare the pattern complex- 
ity for the real and modified topologies. Note that obvi- 
ously, our dynamics has no immediate connection to real 
metabolic dynamics, it rather serves as a dynamic probe 
on complex networks. We find that - in spite of the 
structural and dynamical abstractions - real metabolic 
network topologies exhibit the capacity to maximally reg- 
ularize an imposed complex dynamics compared random- 
ized topologies with identical degree sequences. 



II. MATERIALS AND METHODS 

We use the metabolic networks compiled by Ma and 
Zeng (MZ) 01 . Their data rely on the KEGG database 
f25l | and incorporate enzyme-catalyzed reactions based 
on genomic analyses and biochemical literature. In this 
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FIG. 1: (a) A graph representation of the largest connected component of the unipartite, substrate-centric metabolic network of 
yeast with N = 448 nodes and L — 564 links, (b) The spatio-temporal pattern for the Q(ko) dynamics on the yeast metabolic 
network after a transient of N time steps has been dropped. For visual clarity, we only show 200 nodes for 200 time-steps, (c) 
Entropy signature of the yeast metabolic network (black) and randomized (red), hierarchized (orange) and anti-hierarchized 
(blue) counterparts [see |23|. for details on these randomization procedures] with degree correlations of 0.05, 0, +0.3, and —0.3, 
respectively. Each point represents the entropy signature of the respective network for a randomly selected initial condition. 



unipartite, substrate-centric representation, metabolites 
(nodes) are connected by biochemical reactions (links) 
whenever the catalyzing enzyme is encoded in the respec- 
tive genome. Compared to other metabolic network rep- 
resentations the MZ networks are relatively sparse. This 
is due to the exclusion of frequently occuring metabolites 
like ATP or NADH from those reactions, where they act 
as current metabolites (or carrier metabolites; cf. Q) 
only [see HH for a detailed discussion of current metabo- 
lites in metabolic networks]. As noted in [24|, the ex- 
clusion of current metabolites is reasonable in terms of 
metabolic pathways - otherwise, for example, the path 
length (number of reaction steps) from glucose to pyru- 
vate in the glycolysis pathway would reduce from nine 
to two (24[. Since we are less interested in mass flow 
and catalytic regulation, and more interested in the in- 
formation processing capabilities of such pathways, the 
substrate-centric representation of biochemically mean- 
ingful pathways of the MZ networks is an appropriate 
choice. However, we performed parts of our analysis on 
other network representations (sec Results) with similar 
findings. From the available data we only use the largest 
connected component of the networks similarly to previ- 
ous studies [3,11, Hi!- Hence, we limit our attention to 
the most prominent part of the metabolic networks (on 
average, the largest connected component comprises 54% 
of the nodes and 66% of the links of the complete net- 
work for the 107 species in the MZ database) and avoid 
dynamical artifacts due to isolated residual subgraphs of 
different size. 

Cellular automata (CA) are a well-known tool from 
complexity theory. Under rather general assumptions 
symmetry arguments reduce the rule space for binary cel- 
lular automata on graphs to a set of few outer-totalistic 
[27] ] automata. The effect of these assumptions (like 
isotropy, locality and linearity) on the range of possible 
CA rules will be discussed elsewhere [111)- Most of the 
rules lead to trivial spatio-temporal patterns, e.g. oscil- 
lations or a steady state like the frequently used majority 



rule [29|, [30|. Here, we choose one of those rules, which 
leads to non steady-state behavior on metabolic network 
topologies. Each node acts as a threshold device: the 
state Xi of node i changes from to 1 and vice versa, as 
soon as the density of l's in its neighborhood exceeds k. 
fl can be formalized to 

(x l (t) , Pi < K 

fi(«) : Xi(t + 1) = \ (1) 
[ 1 - Xi (t) , Pi> K . 

The local density pi can be expressed by pi = 
jj-^2,j AijXj, where di is the number of neighbors (the 
degree) of node i and A denotes the adjacency matrix of 
the graph. The rule has been used previously to as- 
sess information processing capacities of synthetic graphs 
[3TI ] . Complex dynamics on MZ networks are achieved for 
a Ko around 0.3. The exact size of the feasible interval 
depends on the maximal degree in the network d max , 
via l/d max . Our implemented complex dynamics can be 
considered as the continuous processing of perturbations 
by the network. Following the argument of 32| that per- 
turbations can travel in both directions of an irreversible 
reaction, we use undirected graphs as metabolic network 
representations . 

Information theory provides tools to analyze spatio- 
temporal patterns, e.g., the Hamming distance, the mu- 
tual information or the Renyi entropy |33| . Here we apply 
the Shannon entropy [34| and the word entropy [3l| , since 
their combination was found to perform well in separat- 
ing different dynamic regimes [311 ]. We thus characterize 
a graph's capacity to process binary information with the 
entropy signature (S,W) on the basis of the specific up- 
date scheme Q. Entropy values are calculated by analyz- 
ing N time-steps for each of the N nodes after a transient 
time of 9N time-steps. 

In our analysis, we take the mean Shannon entropy S of 
the N individual entropies Si, calculated for each node's 
time series separately, as a measure for the structure of 
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the overall pattern, 



N N 



(2) 



The probabilities p° and p\ denote the ratios of O's and 
l's in the time series of node i. Constant nodes yield 
Si = 0, while nodes with a homogeneous distribution of 
O's and l's, that is ever flipping and irregular nodes, con- 
tribute maximally to S with Si « 1. In the dynamical 
regimes regarded in this study, the value of 5* is a measure 
for the homogeneity of dynamical behaviors on the level 
of single states: Large values of S emerge for an overall 
oscillatory or complex dynamics, while smaller values of 
S indicate the existence of constant time series within 
the patterns. 

The word entropy W serves as a simple and easily appli- 
cable complexity measure of the emerging patterns be- 
yond single time steps. To quantify the irregularity of 
the time series of a single node, we count the number of 
constant words, i.e. time blocks of constant cell states, 
of length I. The probability p\ is the number of words 
of length I divided by the number of all constant blocks 
in the time series of node i. The maximal possible word 
length is simply the length t of the time series analyzed. 
The word entropy W of a pattern is the average over the 
individual time series' entropies Wi and is defined as 



W 



N N / t \ 

i=l i=l \ 1=1 I 



(3) 



Patterns with solely oscillatory or stationary behaviors 
result in W=0 while an overall irregular or complex be- 
havior results in high W. 

We use the term entropy for our observables due to 
their formal definition and due to the application of sim- 
ilar concepts in information theory (34j and the theory of 
cellular automata (the word entropy serves as a feasible 
simplification of the 'block entropy', introduced in (35|). 
Note, however, that the two entropy measures defined 
in eqs. and may not be interpreted in the stan- 
dard thermodynamical way, since we average over the 
individual entropies of an ensemble of coupled dynami- 
cal entities. Note also, that the Wi are not bound to the 
interval [0, 1]. 



III. RESULTS 

The metabolic network of the yeast, S. cerevisiae, in 
the MZ database comprises N=752 nodes and L=777 
links. 7V=448 nodes and L=564 links remain in the 
largest connected component of this network. It is char- 
acterized by a diverse mixture of linear chains of nodes 
and hubs connecting chains and single elements (see 
Fig. QJ,). This topological diversity is reproduced in the 
patterns which emerge from running ^I(kq) on the net- 
work (see Fig. [T})): Constant and oscillatory time series 
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FIG. 2: Average entropy signature of the metabolic net- 
works of 22 species and randomized counterparts. The region, 
where the metabolic network topologies (black boxes) reside, 
is clearly separated from the entropy signatures of their ran- 
domized counterparts (red boxes) of same size, connectivity 
and degree sequence. L randomization steps have been per- 
formed on every network. All shown networks have a link to 
node ratio L/N > 1.22. Error bars indicate the spreading of 
the entropy values due to different random initial conditions. 
Additional investigations showed that the separation in the 
plane is robust against variations of position and length of 
the analyzed pattern. A path between the entropy signature 
of the yeast and its randomized counterpart (dotted line) is 
obtained by stepwise randomization. The species abbrevia- 
tions refer to the identifiers used in the MZ database. 



coexist with irregular time series, where nodes switch 
between alternating and constant behavior in an un- 
predictable manner. A randomized network, where the 
topology is changed but the degree sequence of the graph 
is preserved by always switching the links of two pairs of 
nodes as proposed in [36| , shows patterns which are dom- 
inated by irregularities. These differences in the dynamic 
response manifest themselves in different entropy signa- 
tures of the yeast network and randomized counterparts 
(see Fig. [TJj). The mean entropy differences for fi are: 



AS = S 



rand 



' yeast 



0.089, AW = W, 



rand 



yeast 



0.50. Apparently, the yeast topology is capable of reduc- 
ing the resulting entropies, that is, of regularizing the dy- 
namics imposed on it. We checked that any rule among 
the cellular automata on graphs set of rules [28[ lead- 
ing to complex patterns on the yeast topology confirms 
the enhanced regularizing capacity of the real metabolic 
network. Note that here and in the following all graph 
modifications do not only conserve the overall degree dis- 
tribution but also the degree sequence, that is, the indi- 
vidual degrees of all nodes. Moreover, we ensure that the 
randomized graph is connected after every randomization 
step. 

We probe the networks of all species in the MZ 
database and randomized counterparts with the f2( K o) 
dynamics. We find reduced entropy signatures for the 
real topologies for 99 of the 107 species. The net- 
works of the remaining 8 species are extremely sparse 
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FIG. 3: Effect of minimal topological perturbations in yeast. 
Method: (a) For different random initial conditions, we 
yield entropy signatures for the original unperturbed network 
(black diamonds) and a minimally modified topology (gray 
dots); (b) The shifts for each initial condition (gray arrows) 
are averaged (black arrow); (c) Averaging over the ensemble 
of similarly modified topologies leads to the average entropy 
shift for this type of topological perturbation (red arrow). 
Result: (d) Systematic increase of the entropy shift with a 
single randomization step. Different protocols are performed, 
namely arbitrary randomization, randomization under con- 
servation of clustering Coefficient and average path length 
respectively, hierarchization and anti-hierarchization, inter- 
modular and intra-modular randomization. 



and therefore difficult to properly randomize. A reduced 
(S, W) is still observed, if we exclude constant or triv- 
ially oscillating nodes from the analysis. Not too sparse 
metabolic networks (we choose here a link to node ratio 
L/N > 1.22) can be treated in a statistically robust man- 
ner and cluster in the entropy plane at distinctly smaller 
values than their randomized counterparts. 
The modification of a network in single randomization 
steps but at fixed initial conditions leads to a path within 
the entropy plane. Consequently, this assessment of dy- 
namic performance is highly systematic: Small changes 
in graph topology lead to small changes in the entropy 
signature. In Fig. [2] such a randomization procedure 
connects the metabolic network of the yeast and its ran- 
domized counterpart. The opposite direction, from the 
randomized counterpart towards the region of minimal 
entropies can in principle be followed by an appropriate 
randomization process. In this sense, the (S,W) entropy 
plane - or an adequate set of other dynamic observables 
- may provide an orientation for simulated evolution or 
purposive topological design. 

The randomization process alters a number of promi- 



nent topological properties: The average path length 
and the clustering coefficient [13] are reduced, degree 
correlations [23| are diminished and modular substruc- 
tures [U, H[ are disintegrated. In order to investigate 
whether the optimization of the metabolic networks goes 
beyond this set of topological observables, we perform 
different minimal topological perturbations and assess 
the corresponding entropy shift. 

We determine the entropy signatures of the yeast 
network and a modified graph, where only a single 
modification step has been performed, for different 
random initial conditions (Fig. [3£i). Averaging over 
the initial conditions yields the average entropy shift 
for one specific topological perturbation (Fig. [3Jd). 
Averaging over many topological perturbations results 
in a vector which indicates the mean shift in the entropy 
plane due to a specific modification protocol (Fig. \3p). 
We observe an average entropy increase for a single 
randomization, hierarchization and anti-hierarchization 
step (Fig. [3U), where the latter two modifications, 
when iterated, lead to positive and negative degree 
correlations, respectively [see [23l for an explanation of 
the corresponding algorithms]. In order to retain the 
modular structure of the network we identify modules 
with a path length algorithm similar to the one described 
in [38| and randomize solely links within or between 
single modules. Both modifications reveal an average 
entropy increase compared to the original metabolic 
networks (Fig. [5JI). In the case of intra-modular ran- 
domization, the similarity to a random flip suggests a 
topological optimization within the modules. The small 
shift due to inter-modular randomization indicates an 
optimized structure on the highest modular scale. We 
confirmed this result for another module identification 
algorithm based on the topological overlap of nodes [||. 
Moreover, we confirmed an average entropy increase 
for protocols where a link flip is only allowed if the 
average clustering coefficient is conserved or increased 
due to a single randomization step (Fig. [5Ji). The same 
is true for a randomization protocol which conserves the 
average path length (Fig. [3J1). Again, we ensured that 
the randomized networks remain connected through 
the application of the various randomization protocols. 
For all conditional randomization protocols, a average 
entropy increase can result from a limited number of 
sampled topologies. We checked, that this is not the 
case. Furthermore, the entropy shift increases with the 
number of randomization steps, i.e. the modification 
depth, performed on the topology (data not shown). 
Strongly hierarchized and anti-hierarchized networks 
(with considerably positive and negative degree corre- 
lations of +0.3 and —0.3 respectively), cluster in the 
entropy plane far apart from the original networks and 
slightly below and above randomized networks (see 

Fig. Eh). 

Summarizing, the different topological perturbations 
clearly show that the pronounced regularizing capacity of 
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FIG. 4: Average entropy signature for the 22 species of Fig. [5] 
together with synthetic architectures with the identical size 
and connectivity, respectively. We display data for the origi- 
nal metabolic networks (black), scale- free graphs (pink), ran- 
domized scale- free graphs (blue), small- world graphs (green) 
and randomized small- world graphs (yellow). The small- 
world graphs were generated by rewiring 5% of the links of a 
regular chain. The randomized synthetic networks where ob- 
tained by flipping L pairs of links. Notwithstanding the differ- 
ent network sizes and connectivities, the topological variants 
cluster in highly specific regions in the entropy plane. We av- 
eraged over 100 initial conditions for 100 different synthetic 
network realizations each. 



real metabolic networks is not trivially associated with a 
single topological property within the set of graphs with a 
given degree distribution. Moreover, this property goes 
beyond simple degree correlations and modularity and 
are furthermore persistent on the level of individual mod- 
ules. One can ask if the reduced entropy signature of 
real metabolic networks compared to randomized coun- 
terparts is an artifact of the specific network represen- 
tation we discuss throughout this article. However, we 
verified the regularizing capacity of metabolic network 
topologies for another uni-partite substrate-centric rep- 
resentation, where current metabolites have not been re- 
moved from the system. Due to the existence of these 
highly connected hubs, the 43 networks discussed in [26[ 
arc highly connected and show extremely short average 
path lengths. All systems with sufficient network size 
(N > 200) display a pronounced regularizing capacity 
compared to randomized counterparts and an average 
entropy signature increase for singe randomization steps. 
This is also true for two other networks we explicitly 
checked: a directed version of the yeast network derived 
from the MZ database and an undirected enzyme-centric 
network based on the yeast stoichiometric matrix [39j . 
where two enzymes are connected by a link if the prod- 
uct of the reaction catalyzed by the first enzyme serves 
as an educt for the second. Since a regularizing capacity 
is found for other representation of metabolic networks 
as well, we think that this finding a generic property of 
the network architecture of metabolic processes. 

Passing from a real metabolic network to a randomized 



network samples a subset of graphs with a particular de- 
gree sequence. Another level of comparison is achieved if 
we regard the entropy signatures of simple model graphs 
of the same network size and connectivity, that is, the 
same number of nodes and links: Scale-free graphs ob- 
tained with a BA scheme (lo| analog (where the num- 
ber of attached links to newly introduced node varies to 
achieve the desired number of links), and regular graphs 
(where again the number of next-nearest neighbors differs 
between two and three) after a few rewiring steps (small- 
world graphs) according to the original Watts-Strogatz 
scheme |37|. Fig. U shows the corresponding mean en- 
tropy signature and standard deviation for the 22 species 
from Fig. Obviously, the synthetic model graphs re- 
act in a highly specific manner to the imposed dynam- 
ics: They separate from the metabolic networks and from 
each other and cluster tightly in the entropy plane. The 
scale-free graphs exhibit a slightly reduced word entropy 
compared to the real metabolic network, but an (on aver- 
age) higher Shannon entropy. For the small world graphs 
again both entropies are elevated compared to the real 
metabolic networks. Another important feature of Fig.0] 
is that for the different network types the two entropies 
respond differently How do those model graphs react to 
randomization? For these synthetic networks the entropy 
shift measures the bias of the dynamic response, con- 
tained in the particular construction algorithm. For the 
scale-free graphs, we find a small reduction of the average 
word entropy and Shannon entropy under randomization, 
ff we randomize small- world graphs, and thereby break 
up the inherent regular neighborhoods, the randomized 
networks exhibit a reduced average word entropy and a 
elevated average Shannon entropy [3l|. Thus, an average 
entropy increase, as observed for the metabolic networks 
discussed in this paper, is not a trivial effect of the ran- 
domization procedure. On the contrary, for chains of 
regularly connected nodes, a topological perturbation in 
general leads to a decreasing dynamic complexity. In a 
previous study [4p| . we showed this behaviour for syn- 
chronous and asynchronous update schemes with differ- 
ent observables. Moreover we showed, that both topo- 
logical perturbation (i.e., the rewiring of links) and dy- 
namic noise increase the regularizing capacity of an ini- 
tially undisturbed cellular automaton. 



IV. DISCUSSION 

Characterizing complex networks with dynamic probes 
is a novel and promising approach. The recent progress 
in understanding genetic networks by mapping out the 
information flow through the corresponding graphs con- 
stitutes an excellent example of this line of thought [4l| . 
Similarly, our analysis of metabolic networks condenses 
a variety of topological features into the entropy signa- 
ture of the network, while the topological impact enter- 
ing this quantity is selected with respect to its dynamic 
effect. The method discussed here may also serve as an 
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unbiased probe to assess the dynamic performance of dif- 
ferent network types, like protein-interaction networks or 
transcriptional networks, as well as other natural or tech- 
nical networks. Moreover, using system-specific tailored 
dynamics could enhance the structure shaping process in 
engineering artificial network topologies. 
Our key result is that metabolic network topologies em- 
body the capacity to reliably regularize an imposed dy- 
namics - expressed by a systematically reduced entropy 
signature compared to randomized networks with iden- 
tical degree sequences. The large-scale architecture of 
metabolic networks is designed in such a way that com- 
plex and chaotic dynamics are systematically dampened 
out. We believe that the architectural disposition for 
steady-state dynamics in metabolism manifests in the 
observed regularizing capacity of the metabolic network 



topologies. 

In this investigation we focused on the entropies 
averaged over the whole network. Particularly for a 
direct comparison with standard analysis techniques 
of metabolic networks, we believe that the entropy 
distribution over the network could also be extremely 
informative. The individual responses of the system's 
entities to a dynamic probe might help characterize 
different contributions to the functional state of the 
system, or, in terms of metabolic networks, it might help 
understand the principles governing the reorganization of 
steady-state fluxes under topological perturbations [42j . 

We acknowledge the comments of unknown referees, 
which substantially improved the manuscript. 
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